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Abstract. We study a force-based hybrid method that couples atomistic 
models with nonlinear Cauchy-Born elasticity models. We show that the pro- 
posed scheme converges quadratically to the solution of the atomistic model, 
as the ratio between lattice parameter and the characteristic length scale of the 
deformation tends to zero. Convergence is established for general short-ranged 
atomistic potential and for simple lattices in three dimension. The convergence 
is based on consistency and stability analysis. General tools are developed in 
the framework of pseudo-difference operators for stability analysis in arbitrary 
dimension of the multiscale atomistic and continuum coupling methods. 



1. Introduction 

Multiscale methods for mechanical deformation of materials have been investi- 
gated intensely in recent years. The main spirit of these methods is to use atomistic 
models for regions containing defects, and continuum models in regions where the 
material is smoothly deformed. We refer to the recent review [29] for various meth- 
ods and the book [20] for general discussion of multiscale modeling. 

There are two different ways of coupling atomistic and continuum models. One is 
based on energy, and the other is based on force. The energy-based method defines 
an energy which is a mixture of atomistic energy and continuum elasticity energy. 
The energy functional is then minimized to obtain the solution. The force-based 
method works instead at the level of force balance equations. The forces derived 
from atomistic and continuum models are coupled together. The force balance 
equations are solved to obtain the deformed state of the system. 

From a numerical analysis point of view, one of the key issues for these multiscale 
methods is the consistency and stability of the coupled schemes. Taking one of 
the most successful multiscale methods, the quasicontinuum method [26, 38] for 
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example, one of the main issues is the so called ghost force problem [35] , which are 
the artificial non-zero forces that the atoms experience at their equilibrium state. 
In the language of numerical analysis, it means that the scheme lacks consistency 
at the interface between atomistic and continuum regions [16]. In [30], it was shown 
that the ghost forces may lead to a finite size error of the gradient of the solution. 

The stability analysis for the coupling schemes is so far limited to one dimen- 
sional systems, in which case a direct calculation is possible thanks to the easy one 
dimensional lattice structure and pairwise interaction potential. This is no longer 
the case in two and three dimensions, and the extension is by on means easy. More 
general tools for stability analysis are needed, to address in general the multiscale 
hybrid methods. 

In this work, based on existing ideas in the literature, we formulate a force-based 
hybrid scheme for general short-ranged potentials (with some natural assumptions) 
in three dimension. Wc focus on the numerical analysis of the hybrid method, 
which is a representative of a general class of multiscale methods. The solution 
of the proposed method converges quadratically to the solution of the atomistic 
model as the ratio between lattice parameter and the characteristic length scale of 
the mechanical deformation goes to zero. To the best of our knowledge, this is the 
first convergence result for multiscale methods coupling atomistic and continuum 
models in three dimension. 

The convergence result is based on the analysis of consistency and linear stabil- 
ity. To achieve this, we study the linearized operator in the framework of pseudo- 
difference operators. We obtained the stability estimate combining regularity es- 
timate of pseudo-difference operators, consistency of the linearized operator, and 
stability of the continuous problem. These tools developed will help understanding 
multiscale methods in general. 

Before we present the formulation of the method and the main theorem in Sec- 
tion 1.3, we need some preliminaries and notations. 

1.1. Lattice function £ind norms. We will consider only Bravais lattices (see for 
example [3]) in this work, denoted as L. Let d be the dimension. Let {aj} C R'', 
j = 1, - ■ ■ ,dhe basis vectors of the lattice L, hence 




Let {bj} CR^ j = !,■■■ , 



d be the reciprocal basis vectors, given by 



aj ■ bk 



2ndjk- 



The reciprocal lattice L* is then 
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Denote the unit cells of L and L* as F and T* respectively. 

r = {x € M'' I a; = ^ CjUj, < Cj < 1, j = 1, ■ ■ ■ , d}; 

3 

T* = {x£W^\x = Y^ Cjbj, -1/2 < Cj < 1/2, i = 1, ■ • • , d}. 

3 

For e = l/n,n G Z+, we will consider lattice system eh inside domain VL = 
r c M"^, denoted as fig = Q n eL. Note that the lattice constant is £, so that the 
number of points in fi^ is l/e*^. We will restrict to periodic boundary conditions in 
this work, general boundary conditions will be leaved for future publications. For 
a lattice function u defined on eL, we say it is fie -periodic if 

u{x) = u{x'), M x,x' G eh, X — x' = aj for some j G {!, - ■ ■ , d}. 

In particular, an fie-periodic function is determined by its restriction on O^. Func- 
tions defined on fl^ can be easily extended to fie-periodic functions defined on 
eL. 

We also define the reciprocal lattice associated with fig. Let L* = L* n (r*/£). 
Define a subset of Z'^ given by 

j 

hence L* is given by 

hi = {x eR'' \ x = Y^ Hjbj, n e Kg}. 

3 

For ^ gT/, the translation operator is defined as 

{T^u){x)=u{x + enjaj), for a; e M''. 
We define the forward and backward discrete gradient operators as 

= e-'{T^ ^ I) and D", = e-\l - T^), 
where s = /Xjaj and / denotes the identity operator. It is easy to see D'^ _g = 

We say a is a multi-index, if a e Z'' and a > 0. We will use the notation 

d 

j=i 

For a multi-index a, the difference operator D" is given by 

d 

When no confusion will occur, we will omit the subscript e in the notations T^, 
D~g and D" for simplicity. 
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We will use various norms for functions defined on the lattice fig. For integer 
fc > 0, define the difference norm 

0<|a|<A: xeQs 

It is clear that ||-||e,fc is a discrete analog of Sobolev norm associated with if'^(O). 
Hence, we denote the corresponding spaces of lattice functions as H^{0,) and L'^{^) 
when fc = 0. We also need the uniform norms on the lattice Cl^, given by 

||'f*||L~ = max|M(a;)|, 



E max|(i5«u)(a;)|. 



0<|a|<fc 



In the above definitions, wc have identified lattice function u with its fig-periodic 
extension to function defined on eh, and hence the differences are well-defined. 
These norms extend to vector- valued functions as usual. 

Define the discrete Fourier transform for lattice functions / as 

(1.1) /(O = £'^(27r)-'^/2 ^ e-^f-/(x), ieh;, 
and its inverse as 

(1.2) /(x) = (27r)'^/2 ^ e-«/(0, x&Q,. 

We need a symbol which plays the same role for difference operators that (^) = 
1 -I- Aq(^) = 1-1- plays for differential operators. For e > 0, ^ e L*, let 



and 

Alio = 1 + AIM = 1 + E ^UO = 1 + E I (^) • 

It is not hard to check for any ^ G L* , it holds 

(1.3) cA^(0 < AUO < A\0. 

where the positive constant c depends on {bj}. 

The norm of lattice function can be rewritten as 



(1.4) 



11/11^,0 = (27r)'^E 1/(^)1'- 
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Indeed, using Poisson summation formula, 

€eLj j6L| xefie x'efie 

= E £''^(27r)-'^/*(x)/(x') E e'^'^"""'' 

= E (2^)-'^'l/(^)l' = (27r)-<^||/L%. 

Moreover, notice that for ^ e L*, we have 

Therefore, discrete Sobolev norms have equivalent representations using discrete 
Fourier transform: 

c\\f\\ik<Y.^'mm\'<c\\f\\i,, 

with positive constant c depending on k and {aj}. 

For k > d/2, we have the following discrete Sobolev imbedding inequality [24, 
Proposition 6]: 

||/IU|o < C||/||e,fe, 

where C depends on k and Cl. 

1.2. Atomistic model and Cauchy-Born rule. In this work, we will restrict 
our attention to classical empirical potentials. For atoms located at {yi,- ■ ■ , yw}, 
the interaction potential energy between the atoms is given by 

V{yi, ■ ■ ■ ,yN), 

where V often takes the form: 

V{yi, • ■ • , yjv) = E ^2(yi/e, yj/e) + E ysivi/s, Vj/s, yk/s) + ■■■ ■ 

Here we have omitted interactions of more than three atoms. 

Different potentials are chosen for different materials. In this paper, we will work 
with general atomistic models, and we will make the following assumptions on the 
potential functions V as in [19]: 

(1) V is translation invariant. 

(2) V is invariant with respect to rigid body motion. 

(3) V is smooth in a neighborhood of the equilibrium state. 

(4) V has finite range and consequently we will consider only interactions that 
involve a finite number of atoms. 
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The first two assumptions are general [8] , while the latter two are specific technical 
assumptions. 

In fact, for simplicity of notation and clarity of presentation, our presentation 

will be limited to potentials that contain only two-body and three-body potentials. 
Actually, we will sometimes only make explicit the three-body terms in the ex- 
pressions for the potential and omit the two-body terms. It is straightforward to 
extend our results to potentials with interactions of more atoms that satisfy the 
above conditions, following the discussion on the three-body terms. By [25], the 
potential function V is a function of atom distances and angles by invariance with 
respect to rigid body motion. Therefore, we may write 

V^2(2/i,%) = ^2(|2/i-%f ), 

Vs{yi,yj,yk) = V3 {\yi \yi - yk\^, (vi -yj,yi -yk)) , 

where (•, •) denotes the inner product over M'^. We write the two-body and three- 
body potentials this way to make the formula in our calculations easier to read. 

We assume that the atoms are located at fig in equilibrium, with x denoting the 
equilibrium position {x G O^). Positions of the atoms under deformation will be 
viewed as a function defined over O^, denote as y{x) = x+u{x). Hence, u : fi^ ^ R'' 
is the displacement of atoms. We extend u as an fi^-periodic function defined on 
eL. Denote the space of atom positions y as 

Xs = {y : sL, ^ M'^ \ y = X + u, u Oe-periodic, ^ }. 

Hence, y G satisfies 

y{x) — y{x') — X — x' , \f x, x' G eL, x — x' — aj for some j G {1, ■ • ■ , d}. 

The atomistic problem is formulated as follows. For given / : fi^ ^ M'', find 
y G Xg such that 

(1.5) 2/ = arg min /at(z), 
where 

^atW = E E ^(-.-) w - E /(^)^(^)' 

where 

Here S is the set of all possible (si,S2) within the range of the potential. By our 
assumptions, S' is a finite set. Note that as remarked above, we only make explicit 
the three body terms in the potential. In /^t, e'' is a normalization factor, so that 
/at is actually the energy of the system per atom. 

The Euler-Lagrange equations for the atomistic problem is then 

(1.6) ^at[y](a;) = f{x), X G Oe, 
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where 

-^at[y](a;) = (2aiF(,„,,)[y](x)Z?+y(x) + a3V(,,,,,)[y](x)D+ y(x) ) 

(si,s2)es 

+ D- (2a2y(,,,,,)[y](a;)Z)+2/(x) + a3y(,,,.,)M(x)Z?+y(x))), 
where for i = 1, 2, 3, we denote 

the partial derivative with respect to the i-th argument of V. 

To introduce the continuum Cauchy-Born (CB) elasticity problem [8,21,22], 
wc fix more notations. For any positive integer k, wc denote by W'^'^ifl'.M.'^) the 
Sobolev space of mappings y: fl ^ M.'^ such that ||2/||iy'".p < oo- I^i particular, 
W^'^in-W^) denotes the Sobolev space of periodic functions whose distributional 
derivatives of order less than k are in the space U'{Q). For any p > d and m > 0, 
we define X as 

X = {y:fl^R'^\y = x + v,vG W"'+'^'P{fl; R^) f) W^^ifl; R^), J v = 0}. 

As in [19], we have the Cauchy-Born elasticity problem as: find y € X such that 
(1.7) y = argmin J(z), 

where the total energy functional / is given by 

I{z) = [ ( Wcb{Vv{x)) - f{x)z{x) ) da;, 
where v{x) = z{x) — x and Cauchy-Born stored energy density Wqb is given by 

T^cb(A) = | ^ 

' {si,s2)es 

where for ^ G M'^^'^, 

W(s^,S2M) = V{\S1+ SiAf, \S2 + S2Af, (Si + SiA, S2 + S2A) ) . 

The range S is the same as that in the atomistic potential. We have used the 
deformed position y instead of the more usual displacement field u as variable in 
(1.7) in order to be parallel with the atomistic problem. 

The Euler-Lagrange equation for the Cauchy-Born elasticity model is then 

(1-8) J^cbIvKx) = fix), 

where 

-^cbMC^;) = - V- ( DAWcBiVvix)) ) , v{x) ^ y{x) - x. 
Here DaWcb{A) denotes differentiation of H^cb(^) with respect to A. 

Since we are primarily interested in the coupling between the atomistic and con- 
tinuum region, we will take the finite element discretization 7^ be a triangulation of 
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f2e with each atom site as an element vertex with element size e. The triangulation 
is chosen so that it is translation invariant. The approximation space Xg is defined 
as 

X, = {y€ Wl'^in-^R") I y\T G Pi(T), VT G %}, 
where Pi (T) is the space of linear functions on the element T. 

1.3. Force-based hybrid method. We are ready to formulate the force-based 

hybrid method. 

We take p : ^ [0, 1] as a smooth standard cutoff function. The atomistic region 
corresponds to the zero level set of g: fia = {a; | g{x) = 0}, and the continuum 
region corresponds to the region that q equals to 1: Vt^ = {x \ q{x) = 1}. The 
region in between is a buffer between the atomistic and continuum regions. 

The force-based hybrid method is given as: find y{x) e such that 

(1.9) J^i,y[y]{x) = {l-e{x))J^^t[y]{x) + g{x)J^,[y]{x) = f{x), x € n„ 

where J> is the force from finite element approximation of Cauchy-Born elasticity 
problem (1.7). Due to the choice of g, in the atomistic region ila, the force acting 
on the atom is just that of atomistic model, while in the continuum region fic, the 
force is calculated from finite element approximation of the Cauchy-Born elasticity. 

The proposed scheme works in dimension < 3 for general short-range interac- 
tion potentials. The main result for this work is the following quadratic convergence 
result for the force-based hybrid method. 

Theorem 1 (Convergence). Under Assumptions A and B, there exist positive 
constants 5 and M, so that for any p > d and f G W^^'P{i}) fl W^'^{i}) with 
||/||vi/i5.p < (5, we have 

(1-10) ||2/hy-2/at||e,2 <Me2. 

Remark. While we do not attempt in this work to optimize the regularity assump- 
tion on /, we note that it is easy to relax the assumption to / G W^'^{Q) with 
p > d following the remarks below in the proof. 

Remark. The sharp stability conditions Assumptions A and B will be given in 
Section 3. These assumptions are quite natural and physical. We refer to Section 3 
and also [19] for more discussions on the stability conditions and its link to physics 
literature. 

The proof of Theorem 1, which will be viewed as a convergence result for (non- 
linear) finite difference schemes, follows the spirit of Strang's work [37]. In short, 
consistency and linear stability implies convergence. The heart of the matter lies 
in the analysis of consistency and stability, which will be the focus of the proof. 

The rest of the paper is organized as follows. In the next subsection, we review 
some related works. Section 2 discusses the consistency of the scheme. The linear 
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stability is proved in Section 4. The stability estimate is based on the regularity 
estimate of finite difference schemes in Section 3, which is established in the frame- 
work of pseudo-difference operators [9,27,40]. With the preparation of consistency 
and linear stability analysis, the proof is concluded in Section 5. 

1.4. Related works. Recently there are a lot of papers discussing various atom- 
istic/continuum coupling strategies as summarized in the recent reviews [10, 15,29, 
33] , we will only mention some of the works that are closely related to ours and 
refer the readers to these reviews and the references therein. 

The hybrid method resembles several methods in the literature. The most closely 
related method is the quasicontinuum (QC) method [26,38], which is among the 
most popular methods for modeling the mechanical deformation of crystalline solids. 
The QC method contains following ingredients: decomposition of the whole domain 
into atomistic and continuum regions, with the defects covered by the atomistic re- 
gions; degree reduction by adaptive selection of representative atoms (rep-atoms), 
with fewer atoms selected in regions with smooth deformation; and the application 
of the Cauchy-Born approximation in the continuum region to reduce the complex- 
ity involved in compiiting the total energy of the system. 

Both the proposed method and QC method couple atomistic models with non- 
linear Cauchy-Born elasticity model. In some sense, the proposed method can be 
viewed as a smoothened modification of the force-based QC method. Indeed, the 
original force-based QC method amounts to take £> to be a characteristic function 
(so that there is no buffer region). The force-based QC is free of ghost force, and 
it was proven in [12,31] that, for one-dimensional problem, the force-based method 
converges quadratically. However, its convergence behavior remains open for high 
dimensional problem. As will be proved later in the paper, the proposed method is 
stable and also converges quadratically in three dimension. For the understanding 
of the original force-based QC, this work may also provide some new tools and 
insights. 

The Arlequin method [5, 7] and the bridging domain method [6] also adopt a 
smooth transition between atomistic and continuum regions. The difference with 
the proposed scheme is however these methods are energy-based, so that the mixing 
is done at the energy level, while the current method is force-based. Moreover, these 
two methods enforce the consistency between the atomistic and continuum regions 
by imposing certain constraints, while there is no such constraints in our method. 
These methods are suffered from ghost force problems as shown in [29], while the 
proposed method is consistent at the interface. 

The proposed method also shares certain common traits with the concurrent 
Ate coupling method (AtC) proposed in [4]. The AtC method also uses a smooth 
transition between atomistic and continuum regions and is force-based. However, 
the proposed method differs from AtC in the following aspects: (1) our method 
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employs Cauchy-Born elasticity while AtC uses linear elasticity and (2) our method 
is free of ghost force while AtC is plagued by ghost force as demonstrated in [29] . 

Most of the analysis of these multiscale methods limits to the quasicontinuum 
method. In [19], the Cauchy-Born rule for crystalline solids is verified under sharp 
stability conditions. In the language of QC, the authors in [19] actually proved 
the convergence of local QC (the whole computational domain is treated as local 
region). Explicit convergence rate for the local QC can be found in [17,18]. 

For the QC method couples together atomistic and continuum models (nonlocal 
QC method in short), the error estimate can be found in [13,30] and the references 
therein. All these works dealt only with the one dimensional problem, and moreover, 
except [30], the analysis was limited to quadratic potential models, so that the 
system is linear. 

To the best of the authors' knowledge, there is no analysis for the nonlocal 
QC method or other coupling schemes for high-dimensional problems with general 
potential (usually, many-body potential function). The main difficulties lie in the 
analysis of the consistency and stability. For one-dimensional problem, the lattice 
structure is very simple and the pairwise potential function can be handled by a 
direct calculation. However, such an approach cannot be easily extended to high- 
dimensional problem with general potential because the lattice structure and the 
potential function for high-dimensional problem is much more involved. One of the 
main contributions of the current paper is the development of general tools for the 
analysis of consistency and stability. 

Finally, we remark that in this work the analysis of the proposed method, espe- 
cially the stability analysis, is based on analysis of finite difference schemes. The 
readers might wonder why the analysis is not done in the framework of finite ele- 
ment method, as after all, we are dealing with static problems, the systems to be 
solved are "elliptic"; and moreover, the continuum region is discretized by finite 
element method. The reason actually lies in the atomistic part, since the force 
balance equations derived from energy of discrete lattice systems are intrinsically 
of finite difference type. To the best of our knowledge, there has not been yet a 
successful way to put the atomistic equations into the framework of finite element 
analysis. Therefore, to be consistent, we view the finite element approximation in 
the continuum region also as a finite difference approximation. The proof hence 
relies on the analysis of finite difference schemes. This may give a reminiscence 
of the early history about finite element analysis, during when the finite element 
method was also analyzed in the framework of finite difference schemes [36] . Since 
the theory of adaptive mesh is well-established for finite element method, it is an 
interesting question whether one can adopt the finite element analysis framework 
to analyze these multiscale coupling methods. 
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2. Consistency 

We study the consistency of the force-based hybrid method in this section. The 
key is the following lemma, which is a refined version of [19, Lemma 5.1]. 

Lemma 2.1 (Consistency of Cauchy-Born rule). For any y = x + u{x) with u 
smooth, we have 

(2.1) - ^CB[y]|U~ < Ce^uWw^.,^, 

where the constant C depends on V and \\u\\l^, hut is independent of e. 

Remark. The consistency estimate is presented in the form of (2.1) for later use 
in the proof of Proposition 4.2. A bound involves less order of derivatives of u is 
possible, in fact, it is not hard to see from the proof that we have 

(2.2) ||J-at[y]-^CB[t/]||L~ <C£', 

where C depends on V and ||u||^y6,oo. The price is however the dependence of C 
on llMllvvAe.oo is nonlinear. 

Proof. For any x Gflg, and for i = 1, 2, Taylor expansion at x gives 
DtM^) = VlM^)+eVl[y]ix)+e^R2,sAy]ix), 

where, for convenience, we have introduced the short-hands for the Taylor series 
and its remainder: 

Vi[2/](x) = ^(si-V)Ma;), 

Rk,sAy]{x)= [ {k + m-t)''V^+^y{x + etSi)dt, fc e N, 
Jo 

provided that the terms on the right hand side are well defined. Obviously, we may 
write 

(2.3) £>+ = Vj, + e^l + e^R2,s, , £>" = - 6^% - e^R2,-s, • 
For i = 1, 2, 3 and t G [0, 1], let 

Fi{t) = diV^sus.) {\tDtM^) + (1 - i)(si • W)y{x)\\ 

|^£>+y(ar) + (l-^)(s2•V)y(ar)|^ 
{tD+y{x) + {l-t){sr-V)y{x), 

tDly{x) + {l-t){s2-V)y{x))). 

Using Taylor expansion, we get 

(2.4) F,(l) =F,(0)+^^/(0)+i?i[F,](0). 

Here for : [0, 1] — > M, we have introduced a similar short- hand for the remainder 
Rk[Fm = £ ^^p^V'=+iFi(t) dt. 
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Notice that by definition we have 

= 5i^(si,s2)[y](a;); 

Fi{0) = • V)y{x)\\ \{S2 ■ V)y{x)\\ {{s, ■ V)y{x), (s2 • V)y{x))) 

= 5iW^(.„.,)(V^.(x)). 

Therefore, we can rewrite (2.4) as 

9iW(s,,s^){Vu{x)) + sajdijW^s,,s2)i^u{x)) 
+ {e%d,jW^s,.s,){^u{x)) + R,[Fi]{0)) 

Sl,S2 

2i,(si,S2)[Vw](x), 

where for j = 1, 2, 3, 

aj = 2{isj-V)y,Vl.[y]){l-Sjs) 

+ ( ((51 • V)2/, V^^ [y]) + {{S2- V)y, V^^ [y]) ) Sjs, 

bj = 2({sr^)y,^%[y]){i-Sjs) 

+ ( <(5i • W)y, [y]) + {{S2- V)y, Vf, [y]) ) ^,-3, 

c,=2((s,.V)y,i?2,«,[2/])(l-^i3) 

+ ( ((si • V)y, R2,S2 [y]) + {{S2 ■ V)y, R2,s, [y]) ) <5,-3. 

Substituting the equations (2.3) into J^jt [y] (a;) , we obtain 
(si,s2)e5 

+ a3V^(.,,s,)M(Vi^ + eV^^ + 
+(V^, -eV^, -£'i?2,-.j{2a2y(.„.,)[y](V^, +£V2^ +£2i?2,sJ[y] 

+ 9sV^su,,)[y]{Vl+eVl+e'R2,s2)[y]}. 



9^V(^s^,S2)[y]{x) = 

(2.5) 
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Next substituting (2.5) into the above equation, we have 

^Ay\{x) = 

(si,s2)es 

+ Q3,(si,.2)[Vw](Vi, +£V^^ +£'i?2,.2)[2/]} 
+(V^2 -^V^2 -^'^2,-a2){2Q2,(«i,.2)[Vw](V^, +eVl+e-'R2,s.M 

+ Q3,(si,.2)[Vt^](v^i +£V^^ + e'^2,.J[y]}. 

Collecting the terms of the same order, we get 

(2.6) T^t[y\{x) = Co[u]{x) + sCi[u]{x) + s''C2[u]{x) + 0{e^). 

If we change e to — e, the left-hand side of (2.6) is invariant, then the terms of odd 
power of £ in the right-hand side of (2.6) automatically vanishes. Therefore, we 
have 

J-at[y](x) = £oM(x) + e''C2[u]{x) + 0{e^). 
The explicit form of Cq can be written as 

Co[u]{x) = -2(si • V)[(si + (si • V)M)aiW^(,,,,,)(V«(x))] 

- (si ■ V) [(S2 + (S2 • V)^i)53H^(.,,,,)(V^i(x))] 

- 2(S2 ■ V)[(S2 + (.S2 • V)'«)92W^(.i,.,)(V^i(.T))] 

- (S2 • V) [(si + (si • V)u)d^W(,,^,,){Vu{x))\ . 

We see that Cq is the same as the operator that appears in the Euler-Lagrangian 
equation of (1.7). 

The proof of that £2 is of divergence form is similar. Actually, £2 is a quasilinear 
operator, which actually counts for the linear dependence on ||u||vki6.~ on the right- 
hand side of (2.1). To prove (2.1), it remains to estimate terms of 0(e^), which is 
a combination of terms of the form: for q, /3 = 1, 2, 

VL (9^W^(,„,,)(V«)Vi^«), Z + fc = 4,Z,fcGN, 

VL (a.7■3^JW^(.l,.,)(V^i)Vi^^^) , ? + fc = 3,?,fceN, 

We only give the estimate for the first term, and the other two can be bounded 
similarly. Due to chain rule and to Leibniz's rule, Vg^ ( 9iW^(s^^s2)(^^)^s3^ ) is a 
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linear combination of the form 

i=l ^ ^ 

where 7 € N"^ arc mnltiinclccics with I7I — X]i=il7il ItI — 3- Here 

Pi = |Sl + (Sl•V)u|^ P2 = |S2 + (S2•V)«|^ P3 = (Sl + (Si • V)W, S2 + (S2 • V)U) . 

Using chain rule once again, we get, for i = 1, 2, 3, 

\\{Sa ■ V)Pi||Loo < + ||Vu||ioo)||V2u||ioo, 

||(Sa • V)2Pi||ioo < ( (1 + ||Vu|Uoo)||V^u|Uoo + llV^ullioo ) , 

\\{Sa ■ V)'Pi||Loo < ( (1 + ||Vu|Uoo)||V^u|Uoo + ||V2«||ioo||V'u||ioo ) . 

Using Gagliardo-Nirenberg inequality [32] , 

IIV^'^IU- < C\\V"'u\\'l::\\ut-J''^, 0<j<m, 

we have 

\\{Sa ■ V)'=Pi|Uoo < ( ||w|Uoo||V'=+2«||ioc + WV'^+'uWloo ) . 

Using the above inequality, we conclude 

||T|U. < C max ||a^W^(,,,,,)(V«)|U~||V4-l^l«||^«= 

2<|7|<4 

{3 3 
(1 + ii^iiioo) nil v^^+\iu» + nil v^^+^niuoo 
i=l i=l 
3 

+ (l+l|w|lioo) Yl l|V^^+'w||ioo||VT^+2u||ioo||VT'= + lu|Uoo 

i,j,fc=l 



3 

+ (l+||u||ioo) J2 l|V^*+'w||ioo||VT^+lu||ioo||VT'' + lu|U. 




Invoking Gagliardo-Nirenberg inequality again, we obtain 

IITIIioc <C(||u||ioo + ||w||ioo)||Vl°w|Uoo 

+ C{\\u\\U + \\u\\U)\\W\\\loo 

+ C(||u||ioc + ||u||ioo)||V«u||ioc 

+ (7||w||ioc||V^w|Uoo 

6 
i=3 
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Proceeding along the same line, we can obtain the similar bounds for the higher- 
order terms, while ||m||v^i6,oo arises from the following term 

Ra-so. (SiW^(si,S2)(Vu)-R/3,S3b]) • 

Summing up all terms of 0{e^), we get (2.1). 

□ 

Corollary 2.2 (Consistency of finite element discretization). For any y = x + u{x) 
with u smooth, we have 

W^M --^cb[2/]||z,~ < Ce^\\u\\w^e,^, 

where the constant C depends on V and hut is independent of e. 

Proof. The corollary follows Lemma 2.1 by the observation that wc can view the en- 
ergy functional of the finite element discretization as a particular choice of atomistic 
potential energy. 

To be more concrete, let us consider the case = 2, so that each element T e 7^ 
has three vertices. It is straightforward to extend the argument below to higher 
dimensions, with certain complication of notations. 

Let j/e G be the approximation of y so that ye{x) = y{x) for any x € fl^- Let 
Ue = Ue — X. Obviously, we have Us{x) = u{x) for any x e fig. 

Now, for each T £ Ts, Vm^It is a linear function of y^ on the vertices of T. 
Denote the three vertices of T as xo,xi,X2, and Si = {xi —Xo)/s, S2 = {x2 —Xo)/s, 
then VwsIt is the solution of the linear system 

si+ SiA = D+..^ye{xo), 
S2 + S2A = D+^^ys{xo). 

Therefore, let us denote 

VueIt = A(^,^^,^){y,{xo)/e,y,{xi)/e,y,{x2)/e) 

as the solution of the above system. Notice that due to linearity, the map ^(si,s2) 
is independent of e. Hence, for x € T, we can write 

WcsCVUsix)) = WcB{A,^si,S2){ye{xo)/s,ye{xi)/s,ye{X2)/s)) 

(2.7) 

= W"FE,(si,s2)(y£(^^o)/£,ye(a;i)/e,y£(a;2)/£), 

where VFfe,(si,s2) = ^cb ° ^(si,s2)- Denote Sfe as the set of all pairs (si, S2) such 
that {xo, xq + esi,xo + es2} forms the vertices of an element T gTs containing xq 
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(it is easy to sec that Sfe is independent of e). Then, using (2.7), we have 

/ WcB{yUe{x)) 

Jq 

= ^ \T\Wcb{Vu,\t) 

1 d,rr ITT. f yejx) ye{x + esi) Vejx + es2) \ 

= 3! M(si,a2)l>^FE,(si,S2) I — ^> . ) 

^ d ST' Sr T/ f Vsix) ye{x + esi) y^ix + es2)\ 

= 2^ VFE,{sr,s,){^, , I, 

xefie (si,S2)GSfe 

where yFE,(si,s2) = I^(si,s2)|W^fe,(si,s2) andT(s^ j,^) is the triangle formed by vectors 
Si and S2. This indicates that we can view the energy functional in the finite 
element discretization as a particular atomistic potential model, given by three 
body interactions VpE,{si,s2)' by identifying the value of y on nodes as the deformed 
atom positions. 

It is immediately clear that the Cauchy-Born energy density corresponding to 
the atomic potential constructed is just Wcb- Indeed, for a homogenously deformed 
system with deformation gradient A, by definition, the energy of the system is just 
WcB(^)|f^|, and hence the Cauchy-Born energy density is given again by T4^cb(^)- 

With this viewpoint of the finite element discretization as an atomic potential, 
we obtain the conclusion as an immediate corollary of Lemma 2.1. □ 

CoroUciry 2.3 (Local truncation error). For any y ~ x + u{x) with u smooth, we 
have 

(2.8) ll-^hy[y] - J-atbllkr < Ce-'Mw^^.o., 
and 

(2.9) ||J-hy[y] - -FcbMIIl- < Ce^\\u\\w.e^^, 

where the constant C depends on V and \\u\\l°^, but is independent ofe. 
Proof. The inequality (2.8) follows from Lemma 2.1, Corollary 2.2, and 

my[y]-^.t[y]hr = MxKTMi^) ~ ^AvKmLr 

< ll-FatM -.FcbMIIl.^ + - J-cbMIIl-, 

where we have used q{x) S [0, 1]. Similarly, (2.9) follows from Lemma 2.1, Corol- 
lary 2.2, and 

||.Fhy[y] - J-cbMI|l~ < Mx){j^e[y]{x) - J'cBiyKxmLr 

+ 11(1 - p{x)){J^Mx) - •^CB[y](ar))||L~ 

< llJ^eb] -J^CbMIIl- + ||^at[2/] -.FcbMIIl-. 

□ 



CONVERGENCE OF A FORCE-BASED HYBRID METHOD 



17 



3. Regularity estimate 

To analyze the stability property of the proposed force-based hybrid method, 
we use the framework of pseudo- difference operators [27, 40] . In this section, we 
will establish regularity estimate Theorem 3 for the force-based hybrid method. 
This will be one of the key ingredients used to prove stability estimate in the next 
section. 

We study the linearized operator of J-hy. Let us denote "Hhyiw] the linearization 
of Jriy at state u: 

'Hhy[l' 



Sy 

t 



y—x+u 

SO that y-hy[u\ is a linear operator on lattice functions w, given by 

dThy [x + u + tw] 



t^o dt 

It is convenient to rewrite Hhy in the form of a pseudo-difference operator as 

^hyM = hhy[u]{x,fj.)T^, 

where the coefficient hhy[u]{x, fx) is a d x d (probably asymmetric) matrix for each 
X and j2 € A, given by 

(3.1) {hhy[u])al3{x,l^) " 



y—x-\-u 



d{Ti^y)p{x) 

where a, j3 = 1, - ■ ■ ,d are indices. Here A is range of the pseudo-difference stencil 
(note that e ^), which is finite by assumptions. By the definition of ^hy, we have 

(3.2) Ky[u]{x,fi) = (1 - Q{x))h^t[u]{x,ii) + Q{x)he[u]{x,ii), 

where /latM and he{v\ are given by similar equations as (3.1) by replacing J^hy to 
J^at and Fe respectively. 

Define h\^y[u]{x,^) as the symbol of the pseudo- difference operator 'Hhyfu] given 

as 

h\iy[u]{x,i) = ^ /ihyM(a;,/x)exp(z£^/ijaj- • ^) for ^ e L*, 

and similarly for hs[u\ and /iat[w]. By definition, we have for any a; G O^, 
(WhyMefee'^-«),(ar) = (VM)jfe(a;, Oe"'^ 

for 1 < j,k < d and similarly for he[u] and /latM. Here {e/j} are the canonical basis 
of M''. It is also clear that (3.2) implies 

(3.3) hhy[u]{x,^) = (1 - g{x))Kt[u]{x,^) + e{x)h,[u]{x,^). 

In the case that we linearize around the equilibrium state m = 0, we will simplify 
the notation as 

^hy = Uhy [0] , hhy = hhy [O] , hhy = hhy [O] , 
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and similarly for those defined for atomistic model and finite element discretization. 
We observe that by the translation invariance of the total energy /at at the state 
u = 0, 

hat {X,ll) = Kt {fl) , he{x,IJ.) = he (/U) . 

The coefficients are independent of position x, and hence similarly for /lat and he- 

Wc also dciiiotc? HcB as the linearization of Jx;b at the equilibrium state u = 0, 
and define /zcb = ^CB(a^,'^) as its symbol. Note that due to the periodic boundary 
condition assumed on fi, ^ here only takes value in L*. Again, due to the translation 
invariance of the total energy, hen is independent of x. 

Let us start the analysis with the operator ?^hy First, we show that the matrix 
/ihy is Hermitian. 

Lemma 3.1. The matrices /iat(0' hs{^) and hence hhy{x,^) are Hermitian for 
any £ > 0, a; e and ^ e L* . 

Proof. It suffices to prove the result for /lat(C)) as the argument for he{^) is the 
same and the conclusion for hhy{x,^) follows immediately from (3.3). 
Since {J^a.t[y])a{x) = -dht[y\/dya{x), we have 



(/lat)a/3(Ai) 



^ya{x)^{T^^y)p{x) 



dya{x)dyi3{x + eiijaj) 
d-'htly] 



d{T ^'y)a{x + eiijaj)dy0{x + enjaj) 



d{T-i^y)c,{x)dyp{x) 



= (/lat)/3a(-/i), 



j/=x 



where the last line follows from translational invariance of the unperturbed system. 
Therefore, 



A* 3 

= ^(/lat)/3a(-/^)exp(z£^(-/Xj>j • {-£)) 
1^ 3 

= (j2{hs.t)i5a{-n)exp{ieY^{-fij)aj ■ ^) 



for any ^ e L* , where we have used the fact that h^t are real matrices. This proves 
the Lemma. □ 



We make the following stability assumptions about the atomistic potentials, the 
finite element discretization of the Cauchy-Born elasticity model: 
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Assumption A. /iat(0 is positive definite and tlicrc exists aat > such that for 
any e > and any ^ G L* , 

det /iat(0 > a.tK'!eiO- 

Assumption B. he{^) is positive definite and there exists a > such that for any 

e > and any ^ G L* , 

dct/iao>«Ao'(o- 

The Assumptions A and B will be assumed in the sequel without further indi- 
cation. 

Remark. These assumptions are quite natural and physical. In fact, Assumption A 
is just the phonon stability conditions (for simple Bravais lattice) identified in [19] 
represented using the notions of pseudo-difference operators. Assumption B is the 
usual stability condition of a finite element discretization of continuous problem de- 
rived from the Cauchy-Born rule. We note that as a consequence of these stability 
assumptions, the continuous Cauchy-Born elasticity problem is also elliptic, as in- 
dicated by Corollary 4.3 below. From a mathematical point of view. Assumption A 
and Assumption B can be seen as the uniform cllipticity of the difference operator. 

Next, we prove a lower bound for the symbol /ihy, which is crucial for the regu- 
larity and stability estimates. Let us recall an inequality proved by Ky Fan: 

Theorem 2 (Ky Fan's determinant inequality [23]). Let A, B be positive definite 
matrices, then for any A G [0, 1], 

det{XA + {1- X)B) > (det A)^(detB)i-^. 

Corollary 3.2. For any e > 0, x G fl^ o-nd any ^ G L*, we have 

deiKy{x,Cl > min(a,aat)A^j^(0. 

Proof. This is an immediate corollary of Theorem 2. Since for any x, g{x) G [0, 1], 
we have 

dethhy{x,0 = det((l - gix))Kt{0 + Q{^')hs{0) 

> {detKt{Oy~'^''\deth,iOr^^^ 

> al;'^^^a^^-^Al'M 
>min(a,aat)Agf,(0. 

□ 

With these preparations, we now establish the regularity estimate of the quasi- 
continuum approximation. The regularity of discrete elliptic systems is understood 
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by a fundamental result of finite difference approximation by Bube and Strikw- 
erda [9]. They extended the regularity estimate of Thomee and Westergren [39] 
from single elliptic equation to elliptic systems. 

Let us introduce the regular discrete elliptic system following [9]. The concept 
is parallel to the regular continuous elliptic system [1]. 

Definition 3.3 (Regular discrete elliptic system). For i, j = ,d, let Lij be a 

difference operator with symbol lij{x,£). The system of difference equations 
d 

(3.4) ^LijVj{x) = fi{x), i = l,---,d, 

i=i 

is a regular discrete elliptic system, if there are set of integers {(Ji}f^i and {Tj}^^^ 
such that each L^j is a difference operator of order at most CTj + Tj , and if there are 
positive constants C,^o,so such that 

\deth,{x,0\>CAlP{0 

for < e < £o, C e L*, and maxi<j<d|^j| > ^o, where 2p = Y^ii'^i + Ti)- We will 
say that the system (3.4) is regular elliptic of order (cr, r). 

By Corollary 3.2, we immediately have 

Proposition 3.4. Under Assumptions A and B, the finite difference system 

(3.5) Hi,yv = f 

is a regular discrete elliptic system of order (0, 2). 

For the regular discrete elliptic system (3.5), we have the following regularity 
estimate. 

Theorem 3. Under Assumptions A and B, for any v G H^{fl), we have 

(3.6) ||t;||e,2 < C(||?^hy^;||e,0 + Ht'lle.o). 

The constant C is independent of v and e. 

Remark. Theorem 3 is analogous to the interior regularity estimate for elliptic 
partial differential equations given in [2]. The statement of the theorem is just 
rewriting Theorem 2.1 in [9] using the current notation. We note that in [9], Bube 
and Strikwerda proved interior regularity estimates, which clearly implies the a 
priori estimate for periodic case here. 

4. Stability 

The main theorem we will prove in this section is the following stability estimate. 
Theorem 4 (Stability). Under Assumptions A and B, for any v e H^{fl), we have 

(4.1) \\v\\e,2 < C||^hy^;||e,0. 
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Let US make some remarks about the stability result. In general, wc do not know 
whether a stability estimate like (4.1) is valid for the force-based quasicontinuum 
method in general dimension (see [11, 14] for some study in one dimension). Prom a 
pseudo-difference operator point of view, the continuity in x variable of the symbol 
of the linearized operator is crucial for the validity of the strong stability. This is 
also the main motivation to use a smooth transition function q{x) in the current 
scheme. The strong stability property of the scheme will facilitate the numerical 
solution based on iterative methods. 

We also note that the strong stability is also crucial for the extension of the 
current scheme to the time-dependent case. It plays the role of Garding inequality. 
We will leave this to future publications. 

To obtain the stability estimate from the regularity estimate of Theorem 3, we 
need to eliminate ||w|le.o on the right hand side of (3.6). In spatial dimension one, 
this can be achieved by the discrete maximum principle for the finite difference 
equation. This is however no longer the case for higher dimensions, as then we are 
dealing with an elliptic system. The argument we will use is instead similar in spirit 
to the argument used in [1,34] for passing from regularity estimate to uniqueness 
results for elliptic systems. 

The difficulty however is that a compactness argument as in [34] can not apply to 
the finite difference system, as we need a uniform estimate for different e. Therefore, 
instead of using the compactness, the proof is based on the uniqueness of the 
continuous system from ellipticity, the consistency of the finite difference schemes 
to the continuous system, and the regularity estimate Theorem 3. We note that a 
similar approach was considered by Martin [28]. 

In order to connect the finite difference system with continuous PDE, wc need to 
extend grid functions on fi^ to continuous functions defined in f2. For this purpose, 
let us define an interpolation operator as follows.^ For any lattice function u on 
fie, we define QgU G L^{Cl) as 

(4.2) (Q,u)(a;) = (27r)'^/2 ^ e'-«w(0, a; G O. 

Comparing with (1.2), we know that Q^u agrees with u on Cl^- We have the 
following properties of Qg. 

Lemma 4.1. For k >0, there exists constants Ck, Ck > 0, such that for any u, 

Cft||M||/fj(n) < ||(3eM||ff'=(n) < C'fe||u||//fc(n)- 
Proof. The conclusion follows immediately from definition (4.2) and (1.3). □ 



Usual linear interpolations are not sufHcient for our purpose as we need high regularity of 
the interpolated functions. 
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Let X be a standard nonnegativc cut-off function on R"^, wfiicfi is smootfi and 
compactly supported, with ||x||li = 1- Let Xe be the scaled version 

Xe(x)=£-("'')x(e-"x), 

for some a with < a < 1. The choice of the value of a will be specified later in 
the proof of Proposition 4.4. 

Define a low-pass filter operator for / G L^(il) using Xs as Fourier multiplier: 

LTfiO = {2nr^'mre{0 - (27r)'^/^/(e)x(£"0- 

In real space, convolves / with Xe- Note that, using integration by parts, it is 
easy to see that 

(4.3) \rem<Ck\e"^r'', wkez+, 

(4.4) i2nr^'re{0) = 1. 

Hence, is indeed a low-pass filter. For simplicity of notation, we will denote 

Ue = LeQeUe, 

for lattice function on 0.^. 

We state and prove a consistency result for the linearized operator in terms of 
symbols. 

Proposition 4.2 (Consistency of linearized operator). There exists stq > O'nd 
s > such that for any e < eo and ^, G L*, we have 

\hcB{^,r,)-Ky{^,ri)\<Ce\H + iy. 
Proof. By definition, for 1 < j,k < d, 

{hM^,v)=e'{2n)-''/' ^ e-^«-(V),fc(x,r7) 

where fn{x) — e*^'^ for x G and 

(/icB),fe(C,?7) = (27r)-<^/2 / e-'«-dar(/icB)jfc(?7) 

where we have used in the fact that hcB{x, rj) = hcniv) due to translational symme- 
try. Note that we get from the second line from the first line in the above equation 
using the fact that ^ takes value in L* , so that the integral equals to the sum. 
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Hence, taking difference of the above two equations, we obtain the bound 
\hhy{^,v) - hcB{^,r])\ < C sup WHhyiekfv) - 'HcB(efe/^)||L~. 

l<k<d 

Note that by the definition of linearized operators Hhy and Hcb, we have 

y-hyiekfri) - y-CBiekfn) = ^{I'hyix + t{ekfr,)] - J^cb[x + t{ekfn)])- 

Hence, 

||^hy(efe/^) -'HcB{ekfri)\\Lf = lim ^\\J^hy[x + t{ekfn)] - J^cb[x + t{ekfr,)]\\Lp 

< Ce^ekfjw^e.oo < Ce^ekfjHs < Ce^{l + |^?|)^ 

where s is chosen so that the Sobolev inequality 

||/||wi6.oo(n) < C\\f\\Hs{n) 

holds for any / e H^{fl) (s depends on the dimension). Here, we have used 
Corollary 2.3, noticing that ||tefe/^||i,oo is uniformly bounded for 77 as t — >■ 0. This 
concludes the proof. □ 

The proof of Proposition 4.2 actually gives for any e < eq, x G Cle and 77 G L*, 

(4.5) \hi,y{x,'n) -hcB{v)\ < Ce\l + \ri\y. 

Combined with Corollary 3.2, we get as a corollary 

Corollary 4.3. /icb(^) is positive definite and there exists acB > such that for 
any ^ € L*, 

det/icB(C) > acBAl'^iO- 

Proof. Fixed ^ e L* , take ei sufficiently small, so that for £ < ei , ^ e L* (it suffices 
to take El so small that ei^ G F*). Without loss of generality, we can take ei less 
than So in Proposition 4.2. 

From the continuous dependence of matrix determinants on matrix elements, we 
get from (4.5) that for any £ < £1 sufficiently small, x G 51g 

|detX^,,(x,e) -det/icB(OI < + W- 

Combining the last inequality with Corollary 3.2, we get the desired estimate by 
taking £ — >■ 0. □ 

With these preparations, let us now state the key proposition will be used in the 
proof of Theorem 4. 

Proposition 4.4. For {vs}e>o that G H^{^) md \\ve\\s,2 is uniformly bounded, 

we have 



(4.6) lim WHcbVs - y-hyVeWL^m) = 0. 
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Assume the validity of Proposition 4.4, which we wih come back in the end of 
this section, the proof of Theorem 4 follows a reductio ad absurdum. 

Proof of Theorem 4- Suppose (4.1) does not hold, then there is a sequence of func- 
tions {wk} and Sk > such that 

\\wk\\ek,2^oo, asfc->-oo; 
||'HhyWfe||£;,,,o < c, for all k; 

Wk{x) = 0, for all k. 

Set Vk = Wk/\\wk\\ek,2, we then have 

(4.7) \\vk\\e,,2 = l for all fc; 

(4.8) ||HhyWfe||£,,o ^ 0, asfc-)-oo; 

(4.9) Vk{x) = 0, for ah k. 

Since 



y-CBVk = y-hyVk + {Hc^Vk - 'HhyVk)- 

Since H'^hyVfeHe^.o 0, we have 



\\'HhyVk\\L^{n) ^ 0, as A; ^ 00. 
Moreover, by Proposition 4.4, 



WHcBVk - 'HhyVk\\L^{n) ^0, as fc 



oo. 



Hence ||'HcBVfc||L2(n) — > 0. Note also that the average oivk is zero, since Vk{0) = 
0. By the invertibility of ?^cb on the subspace orthogonal to constant function, 
||^^fe||i,2(Q) -)■ 0, as fc 00, while ||t;fc||efc,2 = 1- It follows then ||ufe||e^,o 0. Indeed, 
since 

INIk,i= E <(OI^1t(OI'<i, 

for any 6 > 0, there exist S > and ki, such that for any k > ki, 

(4.10) Yl \vkm'<S/2. 

On the other hand, due to (4.4), there exists /c2, such that for k > k2 

(4.11) E lk'^fc(efi-k(on<'5/4. 

£eL|^,|€|<H 

Moreover, as ||iJfe||L2 0, there exists ks, such that for A; > fcs, 

(4.12) Y M?)I'<V4. 

«eiL* ,|«|<H 



CONVERGENCE OF A FORCE-BASED HYBRID METHOD 25 

Combined (4.10)-(4.12) together, we have for k > inax(fci, fe, fca), 

ii^^iiLo= E \^^'^\'^^- 

Hence, lim;;_^oo||ffe||£fc,o = 0. From Theorem 3, this imphes 

lim \\vk\U^,2 = 0- 

K—^OO 

The contradiction with the choice of Vk proves the Theorem. □ 

Using perturbation, we may extend the results of Theorem 4 to a deformed state 

u. 



Theorem 5 (Stabihty). Under Assumption A and B, there exists 5 > 0, such that 
for any s > and u, \\u\\y^f2,c^ < S and any v G H^{Q), we have 

(4-13) \\v\\e,2 < C\\Khy[u]v\U,o, 

where the constant depends on 6, but is independent of u, v and e. 

Proof. This theorem follows from a perturbation argument of Theorem 4. Denote 
by vq the solution of 

We immediately have 

-Hhy [0] {v-vo) = { Hhy [0] - -Hhy [u] ) V. 
Using Theorem 4, we have 

\\v - vo\U,2 < C\\{nhy[0] - nhy[u] ) t;||,,o < C||Vu||^yi,oc||t;||e,2. 
By triangular inequality, we have 

\\v\\e,2 < \\vo\\e.2 + \\v - Vo\\e,2 

< C||Hhy[0]«o||e,0 + C|| Vu||,^^i,oo ||z;||,,2 

= c\\ny.y[u]v\u.o + c\\s/u\y.,^\\v\u,2 

<C\\nhy[u]v\\,^o + C'6\\v\U,2, 
which gives (4.13) by choosing 6 = 1/(2C). □ 

We conclude this section with the proof of Proposition 4.4. 

Proof of Proposition 4-4- We work in the Fourier domain. By definition, 

{-HcBVeKx) = E e" «/iCB(rE,Ox(£"0^'i(0- 
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Hence, taking Fourier transform, 

= X] ^Gb{£. - 'n,iDx{£°'ri)ve{'q), 



where 



is the Fourier transform of the symbol with respect to x. 
On the other hand, for the discrete system, we have 



where 



Let us compare the difference between HcBVe and HhyVg. We write 



< 1^1(01 + 1-^2(01, 



where 



It suffices to prove that norms of Ii and I2 both go to zero as e ^ 0. Lcit us 
estimate Ii first. By the smoothness of Xi we have |x(£"0 ~ x{£"v)\ < C'e"!^ — ry|, 
hence 

|/i(OI <C^^" El^-^ll^~'(^)^CB(^-r?,r,)||A2(r;)^;(r;)|. 
Define as 

^(0 = 1^1 sup|A-2(,,)Kcb(0»7)|- 

By the smoothness of hcsix, with respect to x and the fact that Hcb is a second 
order operator, we have |^A~^(r/)/icB(^, J?)! < Cj^l""^"^ uniformly in r]. Hence, 
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e /^(L*) as a function of ^. Therefore, 

.1/2 

< Ce"||(?||;i(L.)||Oel'e|U^(n) 

<Ce"||01|,i(L.)|kellffj(n), 

where the first inequahty results from Young's inequahty. This proves that || 
goes to zero as e — > 0. 

Let us consider I2 next. Take a\ G {a, 1), we break I2 into three parts 

i2(0= ^1(0+^2^2 (0 + ^3(0, 

where 

^21 (0 = l|£|>7re-<'iX(£"C) (f^CBi^ - V,V) -hhy{^ - r],r]))ve{r]), 

r,GLJ 

|r;|>27r£"''i 

^23(0 = l|i|<7r£-«iX(e"0 X] (f^CB{^-V,v) -hhy{^-v,v))ve{v)■ 
|?7|<27r£-"i 

We will control each term: I21 is small due to the decay property of x; -^22 is small 
since ^ and 77 is well separated; I23 is small due to consistency. 

I2i- Define w given by 

We observe that is the Fourier transform of 

W{x) = {HciiQeVe){x) - {Qe{'HhyVs)){x) . 

Hence, ||w;||i,2(n) < C||us||£,2- By (4.3), we have 

for any positive integer k. Therefore, we conclude that ||/2i||l2(q) ^ as 
I22' We have 

(4.14) \h2iO\ <cY \^~'in)hcBi^-v,vm''{v)vM\h-r,\>^e-''. 

+ V)hhyi^ - V, Vml{v)Ve{v)\^\i-^\>^e-''i ■ 
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The argument for the two terms are analogous, and let us focus on the first 
term. Consider (^(^) given by 

^{0= SUp\A-\7l)hcB{^,7j)\. 

Since hcB{x,r]) is smooth with respect to x and "Hob is a second-order 
operator, we have if G Z^(L*) as a function of ^. Hence 

lim||</j(0l|5|>,re-°i||;i(L*) =0. 

Therefore, using Young's inequality, the first term on the right hand side 
of (4.14) is bounded by C||v3(^)l|5|>^£-»i \\ii(i^')\\QEVE\\H^(n), which goes to 
zero as e — >■ 0. Hence, I22 goes to zero in L'^ norm. 
I23- Prom Proposition 4.2, we have 

\hcB{^,7j)-h^y{^,7j)\<Ce^{\r]\ + iy 

for some s > 0. As \t}\ < 27re^"i, we have 

\hcB{^,v)-K{^,v)\<Ce^^-^'^^l 

Therefore, 

Ei^3(0P<c J2 ( E ihcB{^-v,v)-K{^-v,v))vMy 

|£|<7r£-"i |T,|<27r£-"i 
<C Yl E |/^CB(^-»7,»?)-/ihy(^ I' 

< C£4-(2^+rf)«i ^ |^;(,7)|2. 

r;|<27re~°i 

Hence, by choosing ai (and also a) sufficiently small that ai < 4/(2s + d), 
we have ||/23||l2 — >• as e — >■ 0. 

Therefore, to sum up, we have proved both ||/i||L2(f2) and ||/2||L2(f2) go to zero 
as e — >^ 0. The proposition is proved. 

□ 



5. Convergence of the force-based hybrid method 

With the consistency and stability results prepared in the last three sections, we 
are now ready to prove the main result Theorem 1. The proof follows the spirit of 
Strang's convergence proof of nonlinear finite difference schemes [37]. 

As a direct consequence of Corollary 2.3, we have the following 
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Corollary 5.1 (Higher order expansion). Under the same assumptions of The- 
orem 1, there exist positive constants 5 and M, so that for any p > d and f G 
W^^'P{Cl) n with ll/llvFi^'P ^ ^> denote y = x + u{x) with u the solution of 

the Cauchy-Born elasticity problem (1.8), we then have 

\\:F^y[y\-f\\Lr<Me\ 

Remark. Using the remark under Lemma 2.1, the regularity assumption of / can 
be relaxed to W^''p{Q.) with p>d. 

Proof of Theorem 1. We take y be that given by Corollary 5.1. It is easy to see 

/' -H^Aty + (1 - t)y\{x) dt-{y-y)= Tkyb] - J'^yM. 
Jo 

Hence y is the solution of (1.9) if and only if 

/' -Hi^yity + (1 - t)y\{x) dt-{y-y) = f- .Fhy[y]. 
Jo 

For any k G (3/2, 2), we define 

B = {yGX, I \\y-y\\e,2<s^}. 

We define a map T : B ^ B as follows: for any y & B, let T{y) be the solution of 
the linear system 

(5.1) /' n^y[ty + (1 - t)y]{x) dt ■ {T{y) -y) = f- F^y[y\. 
Jo 

We first show that T is well defined. Since 

\\ty + {i-t)y-y\\e,2<A\y-y\k2<e'', 

which gives that for sufficiently small e and d < 3, there holds 

\\ty + (1 - t)y - y\\^2,^ < e'^-"'^ < 5, 

where the constant 5 appears in Theorem 5. It follows from Theorem 5 that the 
problem (5.1) is solvable and 

\\T{y)-y\U,2<C\\f-T^y[y\\\e,o 

(5.2) < C\\f - F^MWefi + C\\TAy\ - •^hy[y]||e,0 

where we have used Corollary 5.1. For sufficiently small e, we have 

||T(y)-y|U,2<£^ 

Therefore, T{y) e B and T is well-defined, which in turn implies T{B) C B for 
sufficiently small e. Now the existence of y follows from the Brouwer fixed point 
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theorem. The solution y is locally unique since the Hessian at y is nondegenerate. 
Let us denote the solution as yhy, we then have from (5.2) that 

(5.3) ||y-yhy||s,2<Ce2. 

Proceeding along the same line that leads to (5.2) and using Lemma 2.1, we get 

(5.4) \\y-yAe,2<Ce''. 

Finally, we conclude that yhy satisfies (1.10) by combining (5.3) and (5.4). □ 
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